function M_sim= Algorithm(param,pm,Data,main)

%M_data = [mean(Data.VL_i);mean(Data.VL_f); Data.Lbar_i;Data.Lbar_f];

%
pm.sigma = param(1);
pm.se = param(2);
pm.theta= param(3:2+pm.I);
pm.bc = param(3+pm.I:2+2*pm.I);
pm.phi = param(4+pm.I:2+3*pm.I);

%sigoptions = optimoptions('fmincon','Display','off');
%s0 = [1;1]; s_lb = [0;0]; s_ub = [2;2];
%sigmastar= fmincon(@(s) (Msig_data./Sim_Sigma(pm,Data,s)-1)'*(Msig_data./Sim_Sigma(pm,Data,s)-1), ...
%                         s0,[],[],[],[],s_lb,s_ub,[],sigoptions);

%pm.sigma = sigmastar(1);
%pm.se = sigmastar(2);
%fprintf('Sigma= %0.2f and %0.2f ',pm.sigma, pm.se);

%Calculate Prices and Wages
ze =logninv(1-(Data.N_f+Data.N_i)./pm.M,0,pm.sigma);
zr =logninv(1-Data.N_f./pm.M,0,pm.sigma);

pm.bp = zeros(size(pm.bc));
for r = 1:pm.I
    xbar = integral(@(x) x.*lognpdf(x,0,pm.sigma),zr(r),Inf) ...
                ./(1-logncdf(zr(r),0,pm.sigma));
    Ee= quadgk(@(e) e.*lognpdf(e,0,pm.se),-Inf,Inf);
    pm.bp(r) = pm.ratio(r).^(1./pm.nu).*pm.bc(r)./(xbar.*Ee).^pm.phi;
end




[p,w,~,E_r,LMat, VMat]= PriceWages(pm,ze,zr);


%Generating the correlation between firm size and fraction of contract
%workers
zgrid = linspace(1,2,1000)'.*zr;
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

if min(E_r)<0 && main==1
        M_sim = 1e10;
%        fprintf('Negative Fixed Costs \n');
else
    
ML_sim = [LMat(:,1)./Data.N_i; LMat(:,2)./Data.N_f];
M_sim = [mean(VMat(1:pm.I));mean(VMat(pm.I+1:2*pm.I)); ML_sim; coeff(1)];

end
%}

%{
pm.sigma = param(1);
pm.se = param(2);
pm.theta= param(3:2+pm.I);
pm.bc = param(3+pm.I:2+2*pm.I);
pm.E_i = param(3+2*pm.I:2+3*pm.I);
pm.E_r = param(3+3*pm.I:2+4*pm.I);
pm.bp = pm.ratio.^(1./pm.nu).*pm.bc;
pm.price = ones(pm.I,1);
pm.wage = 10;

[~,~,~,~,~,LMat,~,NMat] = Counterfactuals(pm,Data,pm.bc,pm.bp,1,Inf);

ML_sim = LMat./NMat;
    
lnze =norminv(1-nansum(NMat,2)./pm.M,0,pm.sigma);
lnzr =norminv(1-NMat(:,2)./pm.M,0,pm.sigma);
Msig_sim = Sim_Sigma(pm,lnze,lnzr);
M_sim = [Msig_sim; ML_sim(:,1); ML_sim(:,2)];

%Loss= (M_data./M_sim-1)' * (M_data./M_sim-1);
Loss= sum(abs((M_sim./M_data-1)));
Loss(isnan(Loss)) = 1e10;
%}
